function Data = SAT_SST_func_load_station_data(P)

    yr_sub_st = P.yr_sub_st;
    yr_sub_ed = P.yr_sub_ed;
    
    % *********************************************************************
    % get a subset of coastal and station dataset
    % *********************************************************************
    dir_list         = SAT_SST_func_IO('GHCN');
    file_load        = [dir_list,'Coastal_station_list_',P.app,'_',P.date,'.mat'];
    load(file_load,'ID_use_station','ID_island_station','ID_coastal_station','dis')
    if P.use_mask == 0
        station_list = unique(ID_use_station,'rows');
    else
        station_list = unique([ID_island_station; ID_coastal_station],'rows');
    end
    clear('dir_list','file_load')

    % *********************************************************************
    % Set directory
    % *********************************************************************
    dir_mon          = SAT_SST_func_IO('GHCN');
    dir_load_station = [dir_mon,'matfiles_',P.app,'_',P.date,'/'];

    % *********************************************************************
    % Load in raw station data
    % *********************************************************************
    yr_st            = 1750;
    date_now         = datevec(now);
    yr_ed            = date_now(1);
    N_yr             = yr_ed - yr_st + 1;
    N_sta            = size(station_list,1);
    clear('is_data','lon','lat','lev','T','station_list_mon','T')
    T                = nan(12,yr_ed - yr_st + 1,size(station_list,1));
    lon              = nan(size(station_list,1),1);
    lat              = nan(size(station_list,1),1);
    lev              = nan(size(station_list,1),1);
    is_island        = nan(size(station_list,1),1);
    is_data          = nan(N_sta,N_yr);
    stations         = nan(N_sta,11);

    for ct_station = 1:size(station_list,1)

        station = station_list(ct_station,:);
        if rem(ct_station,100) == 0
            disp([num2str(ct_station),'. ',station])
        end
        
        try
            clear('data','T_sea_level')
            file_load = [dir_load_station,station,'_mon_V4.mat'];
            data = load(file_load);
            
            % Get metadata ------------------------------------------------
            if data.data_lon < 0
                lon(ct_station)    = data.data_lon + 360;
            else
                lon(ct_station)    = data.data_lon;
            end
            lat(ct_station)        = data.data_lat;
            lev(ct_station)        = data.data_lev;
            stations(ct_station,:) = station;
            is_island(ct_station)  = ismember(station,ID_island_station,'rows');
            
            % Apply correction of elevation -------------------------------
            T_sea_level = data.data_T + 6.8 / 1000 * data.data_lev;
            T(:,data.data_yr - yr_st + 1,ct_station) = T_sea_level;
            is_data(ct_station,:)   = nanmean(T(:,:,ct_station),1);
            
        catch
            disp([station ' has no observations'])
        end
        clear('T_sea_level','station')
    end

    l_use     = ~isnan(lon);
    Data.T         = T(:,[yr_sub_st:yr_sub_ed]-yr_st+1,l_use);
    Data.lon       = lon(l_use);
    Data.lat       = lat(l_use);
    Data.lev       = lev(l_use);
    Data.is_island = is_island(l_use);
    Data.is_data   = is_data(l_use,[yr_sub_st:yr_sub_ed]-yr_st+1);
    Data.stations  = stations(l_use,:);
    Data.dis       = dis(l_use);
    
    % *********************************************************************
    % Assign land fraction
    % *********************************************************************
    dir           = SAT_SST_func_IO('GHCN');
    mask          = ncread([dir,'land_percent2_qd.nc'],'lsm');
    mask          = mask([721:1440 1:720],end:-1:1);
    Data.land_fraction = grd2pnt(Data.lon,Data.lat,mask,.25,.25);

end
